*!Version 4.0 Aug13
program define graphsharecdfsumsto, sortpreserve

*** Program to first calculate QTEs overall,
***     and within group,
***	then get share within each percentile 
***	overall for T and C, 
*** 	graph by group QTEs and also shares for T and C
*** adapted from *!Version 4.0 April 12 f graphsharecdf-4.0.ado
***	Uses version 12 graphics
***
***	Puts graphs in `path' directory
***     Labels legends using `legnd1' through `legnd3'
***
*** This  has the shares for each subgroup sum to one over the while CDF
*** Then, figures out where the thing is constant for the 0 earners, and shows last 5 percentiles
*** of that, along with other measures


syntax varlist(min=2 max=2) [if] [in] [aweight/], /*
	*/ VARS(string) NAMEGraph(string) PATH(string) TITL(string) [LEGND1(string)] [LEGND2(string)] [LEGND3(string)] 

	/* parse varlist */
	local depvar : word 1 of `varlist'
	local treat  : word 2 of `varlist'

	marksample touse 

	keep if `touse'  

	parse "vars", parse( " ") 
	preserve

	/* giving weight var a mnemonic name */
	if trim("`weight'") ~= "" { 
		local weightexp "`exp'" 
		local weight = "[aw =`weightexp' ]" 
		di "weight is `weight'"
  		} /*only want to do this once*/

qui { /* start of quietly loop */

	local numvars: word count `vars'
	di "number of vars is : `numvars'"
	di "weight is: `weight'
	
	* get vars
	local vQ1 : word 1 of `vars'
	local vQ2 : word 2 of `vars'
	if `numvars' == 3 {
		local vQ3 : word 3 of `vars'
		}

	di "`vQ1' name is `legnd1'"

	capture drop centile
	gen centile = _n if _n <=99
	
	* set graph so logo is off while still in stata 10
	graph set ps logo off

	*** overall QTEs
	**** Rerun qte to get cutoffs of centiles 
	**** leave off nojcumul so exactly match usual qte for totals
	qte `depvar' `treat' `weight',  q(1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99) nokd nojcumul
	
	** Rename matrices to save
	mat var17co1 = e(b0)'
	mat var17tr1 = e(b1)'

	*** within subgroup
	if `numvars'==3 {	

		* wshare pscore weighted, unweighted version uses unweighted option

		********************************
		*** Start of control group ***
		********************************
		for any wshare1 wshare2 wshare3: gen X =0 if centile<.
		di "`vQ1' is vQ1"

		** New get overall shares of each category
		sum `vQ1' if `treat'==0 `weight' 
		local cov1 = r(mean) // control overall share of var 1
		
		sum `vQ2' if `treat'==0 `weight' 
		local cov2 = r(mean) // control overall share of var 2

		sum `vQ3' if `treat'==0  `weight' 
		local cov3 = r(mean) // control overall share of var 3

		* reset lastzeroqnt so is missing
		capture drop lastzeroqnt

		forvalues cnt = 1/99 {
			* get previous cut
			capture drop lastcut 

			* first time through
			if `cnt'==1 {
			   	*** starts as 0 earnings obs, so last quantile was zero
				gen lastcut=0
				*** lastzeroqnt is last quantile with zero earnings for group 
				gen lastzeroqnt=.
				}
			else {
				gen lastcut=cut
				replace lastzeroqnt = `cnt' if cut==0
				}
			capture drop cut

			* control group cutoff
			gen cut = el(var17co1,`cnt',1)

			* lastcut==cut, still zeros
			if lastcut== cut & cut==0 {
				*weighted
				su `vQ1' if `depvar'<=cut & `treat'==0 `weight'
				capture drop sh
		 		capture drop ash
			 	gen sh = r(mean)
				egen ash = max(sh)
				*** new divide by control variable 1 mean
			 	replace wshare1 = ash/`cov1' if centile==`cnt'
		
				drop sh ash
				su `vQ2' if `depvar'<=cut & `treat'==0 `weight'
			 	gen sh = r(mean)
				egen ash = max(sh)
				*** new divide by control variable 2 mean
			 	replace wshare2 = ash/`cov2' if centile==`cnt'
				
				drop sh ash
				su `vQ3' if `depvar'<=cut & `treat'==0 `weight'
			 	gen sh = r(mean)
				egen ash = max(sh)
				*** new divide by control variable 3 mean
			 	replace wshare3 = ash/`cov3' if centile==`cnt'

				}

			* non zeros
			if cut!=0 {
				*weighted
				su `vQ1' if `depvar'>lastcut & `depvar'<=cut & `treat'==0 `weight'
				capture drop sh
		 		capture drop ash
			 	gen sh = r(mean)
				egen ash = max(sh)
				*** new divide by control variable 1 mean
			 	replace wshare1 = ash/`cov1' if centile==`cnt'
		
				drop sh ash
				su `vQ2' if `depvar'>lastcut & `depvar'<=cut & `treat'==0 `weight'
			 	gen sh = r(mean)
				egen ash = max(sh)
				*** new divide by control variable 2 mean
			 	replace wshare2 = ash/`cov2' if centile==`cnt'
		
				drop sh ash
				su `vQ3' if `depvar'>lastcut & `depvar'<=cut & `treat'==0 `weight'
			 	gen sh = r(mean)
				egen ash = max(sh)
				*** new divide by control variable 3 mean
			 	replace wshare3 = ash/`cov3' if centile==`cnt'


				}


			}
				

		* graph share of each group 
		label variable centile "Percentile index"
		sort centile

		* summarize
		su wshar* if centile<=98
		su wshar* if centile<=45
		*** make sure lastzeroqnt defined right
		su wshar* if centile<=lastzeroqnt

		sort centile

		** New set new part of legend
		local cov1 : di %3.2f `cov1'
		local cov2 : di %3.2f `cov2'
		local cov3 : di %3.2f `cov3'

	
		di "color `namegraph' 1"
		* do color versions
		** set scheme (s1 color, white )
		set scheme s1color

		* graph of shares for control group, taking out ylabels so is set automatically
		scatter wshare1 wshare2 wshare3 centile if centile >= lastzeroqnt-5 & centile <= 98, sort ///
			yline(1) connect(l l l) lcolor(blue green red black) ///
			plotregion(style(none)) ///
			graphregion(fcolor(white) lstyle(none) ilstyle(none) lcolor(white) ilcolor(white)) ///
			legend(region(lstyle(none)) lcolor(white) label(1 "`legnd1' (mean =`cov1')") label(2 "`legnd2' (mean=`cov2')") label(3 "`legnd3' (mean=`cov3')") rows(2)) ///
			lpattern(l - .) msymbol(i i i) saving(`path'/`namegraph'.gph, replace) title("") b1(" ") ///
			xlab(50 60 70 80 90) ylab(, angle(horizontal))
		graph export `path'/`namegraph'.ps, replace logo(off)
		graph export `path'/`namegraph'.eps, replace logo(off)
		!perl -pi -e 's/w setlinewidth/w 3 mul setlinewidth/' `path'/`namegraph'.ps
		!perl -pi -e 's/w setlinewidth/w 3 mul setlinewidth/' `path'/`namegraph'.eps
		!ps2pdf `path'/`namegraph'.ps `path'/`namegraph'.pdf

		di "black and white `namegraph' 1"
		* do black and white versions
		** set scheme (s1 mono, white no background)
		set scheme s1mono

		* graph of shares for control group, taking out ylabels so is set automatically
		scatter wshare1 wshare2 wshare3 centile if centile >= lastzeroqnt-5 & centile <= 98, sort ///
			yline(1) connect(l l l ) plotregion(style(none)) ///
			graphregion(fcolor(white) lstyle(none) ilstyle(none) lcolor(white) ilcolor(white)) ///
			legend(region(lstyle(none)) lcolor(white) label(1 "`legnd1' (mean=`cov1')") label(2 "`legnd2' (mean=`cov2')") label(3 "`legnd3' (mean=`cov3')") rows(2))  ///
			lpattern(l - .) msymbol(i i i) saving(`path'/`namegraph'-bw.gph, replace) ///
			title("")  b1(" ") xlab(50 60 70 80 90) ylab( , angle(horizontal))
		graph export `path'/`namegraph'-bw.ps, replace logo(off)
		graph export `path'/`namegraph'-bw.eps, replace logo(off)
		!perl -pi -e 's/w setlinewidth/w 3 mul setlinewidth/' `path'/`namegraph'-bw.ps
		!perl -pi -e 's/w setlinewidth/w 3 mul setlinewidth/' `path'/`namegraph'-bw.eps
		!ps2pdf `path'/`namegraph'-bw.ps `path'/`namegraph'-bw.pdf

		********************************
		*** Start of treatment group ***
		********************************
		* reset to zero to do treatment group
		for any wshare1 wshare2 wshare3: replace X =0 if centile<.
		
		** New get overall shares of each category
		sum `vQ1' if `treat'==1 `weight'
		local tov1 = r(mean) // treatment overall share of var 1
		
		sum `vQ2' if `treat'==1 `weight'
		local tov2 = r(mean) // treatment overall share of var 2

		sum `vQ3' if `treat'==1 `weight'
		local tov3 = r(mean) // treatment overall share of var 3



		* reset lastzeroqnt so is missing
		capture drop lastzeroqnt

		forvalues cnt = 1/99 {
			* get previous cut
			capture drop lastcut

			* first time through
			if `cnt'==1 {
			   	*** starts as 0 earnings obs, so last quantile was zero
				gen lastcut=0
				*** lastzeroqnt is last quantile with zero earnings for group 
				gen lastzeroqnt=.
				}
			else {
				gen lastcut=cut
				replace lastzeroqnt = `cnt' if cut==0
				}
			capture drop cut

			* treatmentl group cutoff
			capture drop var17tr1
			gen cut = el(var17tr1,`cnt',1)
			* lastcut==cut, still zeros
			if lastcut== cut & cut==0{
				*weighted
				su `vQ1' if `depvar'<=cut &  `treat'==1 `weight'
				capture drop sh
		 		capture drop ash
			 	gen sh = r(mean)
				egen ash = max(sh)
				*** new divide by treatment variable 1 mean
			 	replace wshare1 = ash/`tov1' if centile==`cnt'
		
				drop sh ash
				su `vQ2' if `depvar'<=cut &  `treat'==1 `weight'
			 	gen sh = r(mean)
				egen ash = max(sh)
				*** new divide by treatment variable 2 mean
			 	replace wshare2 = ash/`tov2' if centile==`cnt'
		
				drop sh ash
				su `vQ3' if `depvar'<=cut &  `treat'==1 `weight'
			 	gen sh = r(mean)
				egen ash = max(sh)
				*** new divide by treatment variable 3 mean
			 	replace wshare3 = ash/`tov3' if centile==`cnt'
		
		
				}
			* non zeros
			if cut!=0 {
				*weighted
				su `vQ1' if `depvar'>lastcut & `depvar'<=cut &  `treat'==1 `weight'
				capture drop sh
		 		capture drop ash
			 	gen sh = r(mean)
				egen ash = max(sh)
				*** new divide by treatment variable 1 mean
			 	replace wshare1 = ash/`tov1' if centile==`cnt'
		
				drop sh ash
				su `vQ2' if `depvar'>lastcut & `depvar'<=cut &  `treat'==1 `weight'
			 	gen sh = r(mean)
				egen ash = max(sh)
				*** new divide by treatment variable 2 mean
			 	replace wshare2 = ash/`tov2' if centile==`cnt'
		
				drop sh ash
				su `vQ3' if `depvar'>lastcut & `depvar'<=cut & `treat'==1 `weight'
			 	gen sh = r(mean)
				egen ash = max(sh)
				*** new divide by treatment variable 3 mean
			 	replace wshare3 = ash/`tov3' if centile==`cnt'
		
				}
			}
				
		
		* graph share of each group
		label variable centile "Percentile index"
		sort centile
	
		* summarize
		su wshar* if centile<=98
		su wshar* if centile<=45
		*** make sure lastzeroqnt defined right
		su wshar* if centile<=lastzeroqnt
				
		sort centile

		** New set new part of legend
		local tov1 : di %3.2f `tov1'
		local tov2 : di %3.2f `tov2'
		local tov3 : di %3.2f `tov3'




		} /* end of loop for 3 vars */

	*** for 2 variables
	else if `numvars'==2 { /* start of loop if 2 variables */

		* wshare pscore weighted, unweighted version uses unweighted option

		********************************
		*** Start of control group ***
		********************************
		for any wshare1 wshare2: gen X =0 if centile<.

		** New get overall shares of each category
		sum `vQ1' if `treat'==0 & `touse' `weight' 
		local cov1 = r(mean) // control overall share of var 1
		
		sum `vQ2' if `treat'==0 & `touse' `weight' 
		local cov2 = r(mean) // control overall share of var 2
	
		* reset lastzeroqnt so is missing
		capture drop lastzeroqnt

		forvalues cnt = 1/99 {
			* get previous cut
			capture drop lastcut

			* first time through
			if `cnt'==1 {
			   	*** starts as 0 earnings obs, so last quantile was zero
				gen lastcut=0
				*** lastzeroqnt is last quantile with zero earnings for group 
				gen lastzeroqnt=.

				}
			else {
				gen lastcut=cut
				replace lastzeroqnt = `cnt' if cut==0
				}
			capture drop cut
			* control group cutoff
			capture drop var17co1
			gen cut = el(var17co1,`cnt',1)
			* lastcut==cut
			if lastcut== cut {
				*weighted
				su `vQ1' if `depvar'<=cut & `treat'==0 `weight'
				capture drop sh
		 		capture drop ash
			 	gen sh = r(mean)
				egen ash = max(sh)
				*** new divide by control variable 1 mean
			 	replace wshare1 = ash/`cov1' if centile==`cnt'

		
				drop sh ash
				su `vQ2' if `depvar'<=cut & `treat'==0 `weight'
			 	gen sh = r(mean)
				egen ash = max(sh)
				*** new divide by control variable 2 mean
			 	replace wshare2 = ash/`cov2' if centile==`cnt'

		
		
				}

			* non zeros
			if cut!=0 {
				*weighted
				su `vQ1' if `depvar'>lastcut & `depvar'<=cut & `treat'==0 `weight'
				capture drop sh
		 		capture drop ash
			 	gen sh = r(mean)
				egen ash = max(sh)
				*** new divide by control variable 1 mean
			 	replace wshare1 = ash/`cov1' if centile==`cnt'

		
				drop sh ash
				su `vQ2' if `depvar'>lastcut & `depvar'<=cut & `treat'==0 `weight'
			 	gen sh = r(mean)
				egen ash = max(sh)
				*** new divide by control variable 2 mean
			 	replace wshare2 = ash/`cov2' if centile==`cnt'

		

				}
			}
		
		
		
		gen wsh21 = wshare2+wshare1
		
		* graph share of each group
		label variable centile "Percentile index"
		sort centile

		* summarize
		su wshare* if centile<=98
		su wshare* if centile<=45
		*** make sure lastzeroqnt defined right
		su wshar* if centile<=lastzeroqnt
		
		** New set new part of legend
		local cov1 : di %3.2f `cov1'
		local cov2 : di %3.2f `cov2'

		di "color `namegraph' 1"
		* do color versions
		** set scheme (s1 color, white )
		set scheme s1color

		* graph of shares for control group, taking out ylabels so is set automatically
		scatter wshare1 wshare2 centile if centile >= lastzeroqnt-5 & centile <= 98, sort yline(1) connect(l l l ) lcolor(blue green) plotregion(style(none)) graphregion(fcolor(white) lstyle(none) ilstyle(none) lcolor(white) ilcolor(white)) legend(region(lstyle(none)) lcolor(white) label(1 "`legnd1' (mean=`cov1')") label(2 "`legnd2' (mean=`cov2')") rows(2)) lpattern(l - .) msymbol(i i i) saving(`path'/`namegraph'.gph, replace) title("")  xlab(50 60 70 80 90) ylab( , angle(horizontal))
		graph export `path'/`namegraph'.ps, replace logo(off)
		graph export `path'/`namegraph'.eps, replace logo(off)
		!perl -pi -e 's/w setlinewidth/w 3 mul setlinewidth/' `path'/`namegraph'.ps
		!perl -pi -e 's/w setlinewidth/w 3 mul setlinewidth/' `path'/`namegraph'.eps
		!ps2pdf `path'/`namegraph'.ps `path'/`namegraph'.pdf
		
		di "black and white `namegraph' 1"
		* do black and white versions
		** set scheme (s1 mono, white no background)
		set scheme s1mono

		* graph of shares for control group, taking out ylabels so is set automatically
		scatter wshare1 wshare2 centile if centile >= lastzeroqnt-5 & centile <= 98, sort yline(1) connect(l l l ) plotregion(style(none)) graphregion(fcolor(white) lstyle(none) ilstyle(none) lcolor(white) ilcolor(white)) legend(region(lstyle(none)) lcolor(white) label(1 "`legnd1' (mean=`cov1')") label(2 "`legnd2' (mean=`cov2')") rows(2)) lpattern(l - .) msymbol(i i i) saving(`path'/`namegraph'-bw.gph, replace) title("")  xlab(50 60 70 80 90) ylab( , angle(horizontal))
		graph export `path'/`namegraph'-bw.ps, replace logo(off)
		graph export `path'/`namegraph'-bw.eps, replace logo(off)
		!perl -pi -e 's/w setlinewidth/w 3 mul setlinewidth/' `path'/`namegraph'-bw.ps
		!perl -pi -e 's/w setlinewidth/w 3 mul setlinewidth/' `path'/`namegraph'-bw.eps
		!ps2pdf `path'/`namegraph'-bw.ps `path'/`namegraph'-bw.pdf

		
	
		
				} /* end of loop for 2 vars */
			} 
/* end of quietly loop */
*	version 12

*	pause on
*	pause

	restore


end
	
